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Abstract 

Using data from more than ten-years of observations with the Akeno Giant Air 
Shower Array (AGASA), we pubhshed a result that the energy spectrum of ultra- 
high energy cosmic rays extends beyond the cutoff energy predicted by Greisen [1], 
and Zatsepin and Kuzmin [2]. In this paper, we reevaluate the energy determination 
method used for AGASA events with respect to the lateral distribution of shower 
particles, their attenuation with zenith angle, shower front structure, delayed par- 
ticles observed far from the core and other factors. The currently assigned energies 
of AGASA events have an accuracy of ±25% in event-reconstruction resolution and 
±18% in systematic errors around lO^'^eV. This systematic uncertainty is indepen- 
dent of primary energy above lO^^eV. Based on the energy spectrum from 10^^'^eV 
to a few times 10^'^eV determined at Akeno, there are surely events above lO^'^eV 
and the energy spectrum extends up to a few times lO^'^eV without a GZK-cutoff. 

Key words: Extensive air showers, ultra high energy cosmic rays, energy 
determination 

PACS: 96.40.Pq, 95.55. Vj, 96.40.De, 95.85.Ry 



1 Introduction 

From ten-years of data collected by the Akeno Giant Air Shower Array 
(AGASA), we have shown that the energy spectrum of primary cosmic rays 
extends up to a few times lO^^eV without the expected GZK cutoff [3]. On the 
other hand, the HiRes collaboration has recently claimed that the GZK cutoff 
may be present with their exposure being similar to AGASA [4]. Ave et al. [5] 
have re-analyzed the Haverah Park events and their energies are reduced by 
about 30% using a new energy conversion formula. Although we have already 
published our statistical and systematic errors in the energy determination in 
related papers [3,6,7,8], it is now quite important to reevaluate uncertainties in 
the energy determination of AGASA events with the accumulated data of ten 
years. The uncertainties due to shower front structure and delayed particles 
far from a shower core are also evaluated and described in some detail. 

The AGASA array is the largest operating surface array, covering an area of 
about lOOkm^ and consisting of 111 surface detectors of 2.2m^ area [9,10]. Each 
surface detector is situated with a nearest-neighbor separation of about 1km 
and the detectors are sequentially connected with pairs of optical fibers. All de- 
tectors are controlled at detector sites with their own CPU and through rapid 
communication with a central computer. In the early stage of our experiment 
AGASA was divided into four sub-arrays called "branches" for topographical 
reasons, and air showers were observed independently in each branch. The 
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data acquisition system of AGASA was improved and the four branches were 
unified into a single detection system in December 1995 [11]. After this im- 
provement the array has operated in a quite stable manner with a duty cycle 
of about 95%, while the duty cycle before unification was 89%. 

In a widely spread surface array like AGASA, the local density of charged 
particles at a specific distance from the shower axis is well established as an 
energy estimator [12] since the local density of the electromagnetic component 
depends weakly on variations in interaction models, fluctuations in shower 
development and primary mass. In the AGASA experiment, we adopt the local 
density at 600m, 5(600), which is determined by fitting a lateral distribution 
function (LDF) of observed particle densities to an empirical formula [7] . This 
empirical formula is found to be valid for EAS with energies up to lO^^eV and 
for zenith angles smaller than 45° [13,14]. The relation for converting S'(600) 
to primary energy has been evaluated so far by Monte Carlo simulations [15] 
up to lO^^eV and is 

E = 2.03 X 10^^ ■ 5o(600) eV , (1) 

where 5*0(600) is the 5(600) value per m^ for a vertically incident shower. 
This conversion relation is derived from electron components for air showers 
observed 900m above sea level. In §3.3, a new conversion constant is evaluated 
taking account of the average altitude. More modern simulation codes have 
been used in this simulation. 

In the southeast corner of AGASA there is the Akeno Ikm^ array [16]. This 
is a densely packed array of detectors covering an area of Ikm^ operated since 
1979. This array was used to determine the energy spectrum between 10^^'^eV 
and 10^^'^eV. In this experiment, the total number of electrons, known as the 
shower size N^, was used as an energy estimator. The relation between this 
energy spectrum and the AGASA energy spectrum is discussed in §4. 



2 Densities measured by scintillation detectors 

The AGASA array consists of plastic scintillators of 2.2m^ area, and the 
light from these scintillators is viewed by a 125mm diameter Hamamatsu 
R1512 photomultiplier tube (PMT) at the bottom of each enclosure box. The 
scintillators are 5cm thick (0.14 radiation lengths). The enclosure box and a 
detector hut are made of steel with 2mm and 0.4mm thickness, respectively. 
There is another type of enclosure box used for 17 detectors in the Akeno 
branch, in which PMT is mounted at the top of each enclosure box. 

In order to cover a dynamic range from 0.3 to a few times 10^ particles per 
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detector, a logarithmic amplifier is used [9]. The number of incident particles 
is determined from the pulse width, which is obtained by presenting such 
a signal that decays exponentially with a time constant of r ~ 10/is to a 
discriminator with a constant threshold level V^. The relation between the 
number of incident particles and the pulse width td is given by 

Vd = Ve-^^l^ = kNe-^^'^ , (2) 



where V is the pulse height {V = kN, k is a constant depending on the gain 
of the amplifier). By defining the pulse width for = 1 as ti, one obtains 

lniV=^^^ . (3) 



Figure 1 shows a typical pulse width distribution (PWD) for omni-directional 
muons, and its peak value is used as ti in the experimental convenience. 

In the Akeno experiment, the original definition of a "single particle" was 
based on the average value {PH^^J of a pulse height distribution (PHD) from 
muons traversing a scintillator vertically [17]. This PH^^^ is accidentally coin- 
cident with the peak value PH^^^j^ of the PHD of omni-directional muons, since 
the PHD is not a Gaussian distribution but is subject to Landau fluctuations 
and coincidental incidence of two or more particles. If we express the pulse 
height corresponding to the peak value of PWD as PW^^^f^, the parameter 
PH%,^ is related to PW'^,,^ by 

PW',,ak = \ [pHlak + {PHIJ^^ . (4) 



This equation can be derived under an assumption that main part of PHD 
is expressed by a Gaussian distribution with a standard deviation of a. By 
converting the pulse height of this Gaussian distribution to the pulse width by 
using Equation (2), and by evaluating ti at the maximum of the distribution. 
Equation (4) is obtained. With PH^^^^ = 1.0 and a = 0.35, PW^^^^ = 1.1. 
The density measured in units of PW^^^j., therefore, is 1.1 times smaller than 
that measured in units of PH^eaki— P^ave)- other hand, the density 

measured with a scintillator in units of PH^y^ times larger than the 

particle density measured with spark chambers between 10m and 100m from 
shower cores [18]. The efficiencies and other details of the spark chamber are 
described in [19]. This factor 1.1 was interpreted as due to the transition effect 
of electromagnetic components in scintillator compared to the electron density 
measured in spark chamber [18]. 

This means that the density in units of -PW£.aA; corresponds to an electron 



4 



density measured by a spark chamber, given that the ratio of densities mea- 
sured with scintillators and spark chambers is 1.1 . The number of particles in 
units of PWp^^^., therefore, coincides with the electron density and has been 
applied so far to estimate primary energy using Equation (1). 

To examine whether a "single particle" is appropriate, the CORSIKA pro- 
gram has been used to simulate densities measured by a scintillator of a 5cm 
thickness [20]. In this simulation, a "single particle" corresponds to PH^^^. 
Figure 5 plots the lateral distribution of energy deposit in the scintillator in 
units of -Pi^Q^e (closed circles), and it is compared with the experimental LDF 
(dashed curve). The simulated LDF is flatter than the experimental one. The 
simulated density reflects the number of electrons near the core (up to about 
200m from the core), but becomes larger than the electron density with in- 
creasing core distance. Recently we have also studied the detector response 
with the GEANT simulation [21,22]. In this simulation, a "single particle" 
is defined as the peak value of log^Q (energy deposit in scintillator) for omni- 
directional muons with their energy spectrum to represent the experimental 
P^peak- Here we take account of the real configuration of a detector, con- 
version of photons in the wall of the enclosure box and the detector hut, 
scattering of particles, decay of unstable particles (pions, kaons and etc), and 
the 4-momentum of shower particles. The shape of the lateral distribution is 
nearly consistent with the experimental LDF, though it is also a little flatter 
than the experimental one. The details will be described in [22]. 



3 Evaluation of uncertainties on energy estimation in AGASA 

3.1 Detector 

The detector positions were measured using a stereo camera from an air- 
plane with accuracies of AX, AY = 0.1m and AZ = 0.3m. The cable lengths 
(the propagation delay times of signals) from the Akeno Observatory to each 
detector is regularly measured with accuracy of 0.1ns in each RUN (about 
twice a day). Figure 2 shows the variation in cable length for a typical de- 
tector as a function of day. A discontinuity around 50,000 MJD is due to 
the movement of the detector position and another one is due to the system 
upgrade in 1995. 

In Equation (3), there are two parameters which should be determined. The 
first one is the "single particle" ti (= PWp^^j.). In the AGASA experiment, 
pulse widths of all incident particles are recorded and their PWD is stored in 
the memory as shown in Figure 1, and then ti is determined in every RUN. 
Figure 3(a) shows the time variation of ti for a typical detector over 11-years 
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of operation. There is a clear seasonal variation with a ±3% fluctuation, but 
this variation has been calibrated in the air shower analysis using monthly 
data. The variance cr'^iti) within each month of data is determined and 
is shown for all detectors in Figure 3(b), with a{ti)/{ti) < 0.7% for a 68% 
C.L. 

The second important parameter is the decay constant r. Although we have 
directly measured the r values with a LED several times during AGASA's 
operation, they are not enough to estimate the time variation of Ar/r. We 
estimate this variation using the observed PWDs. Assuming the density (par- 
ticle number) spectrum of incident particles in a detector is J oc A^""^ , one 
obtains 

Aln/ = -7— — = -7 [td = Kx) . (5) 

iV r 



The ratio a = A In J/ Ax is the slope of the PWD, so that the ratio Ar/r is 
expressed by 

^ = -^ . (6) 
T a 



Figure 4(a) shows the time variation in Aa/a for a typical detector. The 
variance o"^(^) is determined throughout the observation time (11 years) 
for each detector. For all 111 detectors, cr{^) is plotted in Figure 4(b), and 
o'{—)/{—) < 1.6% at a 68% C.L. This fluctuation causes an uncertainty in 
density estimation of 4%, 7%, 11% and 15% for 10, 10^, 10^ and 10^ particles 
per detector, respectively. It should be noted that cr(— ) includes not only 
the change of r but also that of 7 caused by varying atmospheric conditions 
(temperature and pressure). The real variation in Ar/r is, therefore, smaller 
than that for Aa/a discussed above. 



3.2 Air Shower Phenomenology 



(a) Lateral Distribution Function: 

In general, the number of detectors with measurable particle density 
is not enough to treat both the core location and the slope parameter in 
lateral distribution as free parameters for most of the observed showers. In 
order to avoid systematic errors in determination of lateral distribution, 
we have used only showers well within the boundary of the array and 
introduced a new parameter which represents the degree of goodness in 
locating the core position which is estimated by the maximum likelihood 
method. The details is described in [7]. 
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The empirical LDF thus determined [7] is expressed by 



^M-(i^) O^d V +(1000) I • 

where r is the distance from the shower axis in meters. The Mohere unit 
Rm is 91.6m at the Akeno level. The parameter is a function of zenith 
angle 9 expressed by 

r/= (3.97 ±0.13) -(1.79 ±0.62) (sec ^-1) . (8) 

The uncertainty in the energy determination of showers due to the limited 
accuracy in determination of r] was discussed and estimated to be ±10% 
by Yoshida et al. [7]. 

With observed showers, we have confirmed in [14] that the empirical 
formula of Equations (7) can be applied to showers with energies up to 
j^gi9.8gY with core distances up to 3km by improving the errors of 
Equation (8) by 

r/= (3.84 ±0.11) -(2.15 ±0.56) (sec 0-1) . (9) 

We may extend the present lateral distribution up to the highest energy 
observed, since any energy dependence of rj has not been observed. In 
the same manner as Yoshida et al. [7], the systematic effect on 5(600) 
estimation due to uncertainties in Equation (9) is evaluated to be ±7% 
for air showers with zenith angles smaller than 45°. 
(b) Atmospheric Attenuation: 

Since an inclined air shower traverses the atmosphere deeper than a 
vertical shower, a shower density Sg{600) observed at zenith angle 6 must 
be transformed into 5*0(600) corresponding to a vertical shower. The at- 
tenuation of 5*0(600) is formulated as follows: 



Se{600) = 5o(600) exp 



-^(sec0-l)-^(sec0-l)2 
Al A2 



(10) 



where Xq = 920g/cm2, Ai = 500g/cm2 and A2 = 594li^gg/cm2 for 6 < 
45° [7]. The uncertainty in 5(600) determination due to the uncertainty 
in the attenuation curve of 5'(600) was also discussed there. 

The attenuation curve of 5'(600) is now under reevaluation with the 
accumulated data up to zenith angles 9 < 60° and using a modern simu- 
lation code. Equation (10) is valid for observed events with 6 < 45°, and 
AIRES simulation agrees well with the experimental formula and data 
[14]. Another important point of this study is that As are independent 
of energy. The uncertainty in 5*0(600) due to this transformation is es- 
timated to be ±5%; this value is also reduced from Yoshida et al. [7] 
because of the increased amount of observed showers. 
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In the analysis we have assumed the symmetric distribution of particles 
around the shower axis by neglecting the difference of attenuation to 
each detector. The experimental LDF is determined by the average of 
these densities, and the fluctuation around the average is determined with 
these densities. In analysis procedure described in §3.4, this fluctuation 
is included in a, which is used as a weight of each detector. 

(c) Accidental Coincidence: 

Because we use the log-amplifier described above, the density could be 
overestimated if an accidental signal hits on the tail of the exponential 
pulse above the threshold level of the discriminator. The counting rate of 
the scintillation detector (area 2.2m^) is about 500Hz for signals exceeding 
the threshold of 0.3 particles per detector. The accidental coincidence of a 
background particle hitting a detector within the pulse of a single particle 
(10/is width) occurs with a chance probability of 10 x 10~^ x 500 = 5 x 
10-3. In the same way, one obtains 1.65 x 10 ^ for the probability during 
a 10 particle pulse (33yus width) and 5.6 x 10"^ for a 100 particle pulse 
(56/xs width). This means that one of 200, 61, or 18 detectors in each case 
may record larger values than the real density. Frequency of accidental 
pulses is much less than that of delayed particles described in §3. 2(e) and 
their pulse heights are smaller than 1 particle. With the present analysis 
method described in §3.4, a detector which deviates more than 3o" from 
the average LDF is excluded and hence this effect is negligible. 

(d) Shower Front Structure: 

Given our use of the log-amplifier, the density is properly estimated so 
long as the thickness of a shower front is less than a few 100ns. However, if 
the thickness is larger than this time width, we should take this effect into 
account to estimate the incident particle density appropriately. Figure 
6 is an example of the arrival time distribution observed with a 30m^ 
scintillator [6], operated by the Yamanashi university group and triggered 
by AGASA. The core distance is 1,920m and the primary energy is 2 x 
lO^^eV. Not only is the particle arrival time distribution broad, but 5 
particles are delayed more than 3/is in this 30m^ detector. 

The arrival time distribution of shower particles has been measured 
with a scintillation detector of 12m^ area together with the 30m^ detector. 
Signal sequences of arriving particles are recorded in time bins of 50ns for 
the 30m^ detector and 20ns for the 12m^ detector in coincidence with the 
AGASA trigger. The details of these experiments are described in Honda 
et al. [23,24]. From these experiments, the average shape of the arrival 
time distribution is expressed by 



where the scaling parameter is 168ns, 212ns, and 311ns at r = 534m, 
750m, and 1,050m, respectively. These are shown by solid lines in Figure 




(11) 
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7. Beyond 1,050m the events are too few to determine the average to- 
However, if we extrapolate the relation assuming tg oc r, we find to = 
440ns at r = 1,500m and 570ns at 2,000m. If we extrapolate it assuming 
log to oc r, to = 490ns at r = 1,500m and 850ns at 2,000m, respectively. 
The to distribution of each event in these distant ranges seems to be nearer 
to, but shorter than the values extrapolated with log to oc r. We use, 
therefore, the latter values in the present evaluation as upper bound up 
to 2,000m and for events with energies up to lO^^eV. These distributions 
are also drawn with dotted lines in Figure 7. The solid curves in Figure 6 
correspond to the time distribution with to = 800ns, which support the 
extrapolation with log to oc r. 

Using to and the number of incident particles as parameters, we have 
derived the ratio (the overestimation factor) of the estimated density 
due to the broadening of the shower front structure to the density with 
to = as a function of core distance and primary energy. The results are 
drawn in Figure 8. The factor is nearly independent of primary energy 
up to a few 1,000m. The factor increases rapidly with core distance above 
1,500m, but it decreases suddenly again at those core distances where the 
observed number of shower particles is near unity. From this figure, the 
overestimation factor for the density at 600m is +3.5% and that around 
1km is +6% for lO^^eV showers. With our analysis method described in 
§3.4, the present S'(600) may be overestimated by about 5% due to the 
broadening of the shower front structure with its fluctuation about ±5%. 
(e) Delayed Particles: 

As shown in Figure 6, there are particles at large core distances which 
are delayed by more than a few micro seconds with respect to normal 
shower particles. It was shown in the prototype AGASA experiment that 
pulses delayed by more than 4/is are most likely to be low energy neutrons 
with energy 30-40MeV and the fraction of these pulses to the total shower 
particles is a few % between 1km and 3km [25]. Based on this result, we 
have so far assumed that the effect of delayed particles on the ^(GOO) 
determination is within the errors due to other effects. 

In the following, we evaluate the effect of delayed particles on the 
£'(600) determination with accumulated data from ten years of opera- 
tion. If a delayed particle, whose energy loss in a scintillator corresponds 
to Nd particles, hits a detector with time delay tf) with respect to A^'j 
incident particles at t = 0, the pulse height V(t) is expressed by 

V{t) = Ni exp (^-^) + Nd exp (-^i^) = OF ■ N, exp (^-^) (12) 
0..l/fexp(?£) ] 

where Aj and A^^ are in units of a "single particle" , and OF is the over- 
estimation factor due to delayed particles. The OF value depends on 
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Nn/Ni and tr,- It is, therefore, quite important to evaluate the density 
of delayed particles and their energy loss in a 5cm scintillator experimen- 
tally. These values have been measured with the scintillation detectors of 
30m^ and 12m^ described in the previous section. 

Figure 9 shows the ratio of delayed particles (delay time to > 3/is 
and pulse height No > 1.0 particles) to all shower particles measured as 
a function of core distance for three energy ranges: log(Energy [eV]) = 
18.5-19.0 (open circles), 19.0-19.4 (open squares) and above 19.4 (closed 
squares). In the same figure, the previous Akeno result by Teshima et 
al. [25] and the result by Linsley [26] are also plotted by small open 
and closed circles, respectively. In these measurements, ratios of delayed 
particles with delay time > 4/is and pulse height No > 3.0 particles 
to all particles are plotted for showers of energies around lO^^eV. When 
we take account of the different selection conditions it may be concluded 
that the ratio of the number of delayed particles to all shower particles 
depends on core distance, and is almost independent of primary energy 
from lO^^eV to lO^^eV. 

Figure 10 shows the ratio No/Ni as a function of core distance ob- 
served by the 30m^ or 12m^ detectors. Here, No represents the energy 
loss in scintillators of the delayed particles (delay time to > 3/is and 
pulse height No > 1.0 particles) measured in units of PH^^,^, and Ni is 
sum of all shower particles also in units of PH^y^ ^r showers in the same 
distance and energy bin. Open circles represent the ratio for showers of 
energies between 10^^'^eV and 10^^ °eV, and closed squares for energies 
above lO^^-^eV. 

In Figure 11, the delay time and No of a delayed particle (units of 
Pi/°„J are plotted function of core distance. Since there is no ap- 
preciable difference for different primary energies, delayed particles for 
all energy ranges are put together in this analysis. Details of the shower 
front structure and delayed particles will be described elsewhere [27]. 

Using the No/Ni values in Figure 10 and the delay times in Figure 11, 
the OF values are estimated from Equation (12) and plotted in Figure 
12 as a function of core distance. These OF values are independent of 
primary energy. This is understood as follows. For lO^^eV showers, the 
density at a core distance of 1km is 6/m^ (for AGASA detector N^ = 
13.2/2.2m^), so that the density overestimation OF due to delayed par- 
ticles with No = 3 is expected to be 1.31 for to = 3/is and 1.37 for 
to = 5/is. However, the density of the delayed particles is so small that 
only one in 10 detectors around 1km will be hit by delayed particles. 
Since the average density is determined by several detectors, the ^(GOO) 
overestimation is limited to 4%. On the other hand, for lO^'^eV showers 
all detectors around 1km from the air shower core are likely to be hit 
by delayed particles. However, the density overestimation (OF) of each 
detector is 1.04 because the density of shower particles around 1km is 
large (A^^ = 132) at lO^OeV. 



10 



From our analysis procedure described in §3.4, S'(600) may be overes- 
timated due to delayed particles by about +5% ±5%, independent of pri- 
mary energy. It should be noted that the AGASA LDF is consistent with 
that from electromagnetic components and muons simulated as described 
in §3.1. If we include the simulated results on low energy neutrons (de- 
layed particles) using AIRES code, the LDF becomes much flatter than 
the observed LDF beyond 1km from the core. A possible flattening of 
LDF due to delayed particles is not observed experimentally up to 3km 
from the core and up to lO^^eV. 

3.3 Energy Estimator 

The particle density S'o(600) in Equation (1) is evaluated as the electron 
density and the AGASA density in units of PW^^^f^ corresponds to the electron 
density since the ratio of densities measured with scintillators and a spark 
chamber is 1.1 as described in §2. Since this ratio is not measured beyond 
100m, it is necessary to evaluate the conversion factor from 5(600) measured 
in units of the AGASA "single particle" to primary energy. 

The new conversion formula obtained is described in Sakaki et al. [22] and 
hsted in Table 1. In this simulation a "single particle" is defined as PWp^^f. 
in accordance with the experiment. The energy conversion formula of Equa- 
tion (1) was estimated for the 900m altitude of the Akeno Observatory. Since 
the Akeno Observatory is located on a mountain side, core positions of most 
events are lower than this altitude. At the average 667m height of the AGASA 
detectors, the atmospheric depth is 27g/cm^ larger than that at the Akeno 
Observatory. In the new simulation, this altitude is applied. 

If we evaluate the difference in a factor a in Table 1 due to the difference 
of average altitudes 900m and 667m, it leads to a 7% increase at S'(600) = 1. 
That is. Equation (1) evaluated at 900m is revised to 

E = 2.17x 10^^-5o(600)^-° eV , (13) 

at 667m and the result agrees with the factor a calculated using the QGSJET 
interaction model and a proton primary by Sakaki et al. [21] in Table 1. 

In order to see the differences due to simulation codes and hadronic inter- 
action models, the simulation by Nagano et al. [20] using CORSIKA is also 
listed. In this simulation, the density in units of PH^^^j. is used and the av- 
erage altitude is 900m. Taking account of 10% overestimation of a "single 
particle" (Pif^g^^) and the 7% underestimation due to differences in altitude, 
we may directly compare these results with Sakaki et al. [21]. In each simula- 
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tion, the differences are within 10% between QGSJET and SIBYLL hadronic 
interaction models and are within 10% between proton and iron primaries. 
The difference due to the simulation code itself is within 5%. 

It is, therefore, reasonable to use a revised energy conversion formula by 
taking the average of these simulation results at 667m for proton and iron 
primaries with AIRES (QGSJET, SIBYLL), CORSIKA (QGSJET, SIBYLL) 
and COSMOS (QCDJET) yielding 

E = 2.21 X 10^^ ■ 50(600)^-°^ eV . (14) 



That is, the AGASA energies so far published must be shifted by +8.9% at 
2 X lO^^eV, +12.2% at lO^^eV and +13.1% at lO^^eV. The systematics due 
to the simulation codes, interaction models, and mass composition may be 
within 10%. The intrinsic S'(600) fluctuation in shower development is less 
than 6% under each combination of primary mass and interaction model with 
the AIRES simulation [22]. This small difference among interaction models 
and compositions is an advantage of measuring S'(600) using scintillators, in 
which observed particles are dominated by electromagnetic components with 
a small contribution of muons. 

From the above discussion. Equation (1) used so far by the AGASA group 
gives the lowest limit in the conversion from 5(600) to primary energy. It may 
be more reasonable to increase the energies +10% ±12%. A detailed study of 
this topic will be found in [22]. 



3.4 Analysis 



Our analysis procedure for an air shower event is based on an iterative 
process to find the arrival direction of a primary cosmic ray and to search 
for the core location and the local density S'(600). To start, we assume an 
initial core location at the center of gravity of the density distribution of an 
observed event. Next, the arrival direction (zenith angle 6, azimuth angle 0) 
is determined by minimizing the function: 



n 



1 " 

— [{T, - Tfin, e, 0) - T,(p„ R,) - To}' /T,(p„ R,f] , (15) 



where Tj is the observed time of the first particle incident on z-th detector (the 
detector location is expressed by r, measured from the core position and Ri is 
the distance from the shower axis, and pi denotes the observed density). Here, 
Tf is the propagation time of the tangential plane of the shower front, is the 
average time delay of the shower particles from this tangential plane, is the 
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average deviation of shower particles, and Tq denotes the time when the core 
hits the ground. The parameters - T^^siPi, Ri) ^ of shower front structure are 
obtained experimentally [28]. At this step, those detectors that make large 
are excluded in the calculation as signals with accidental muons. This exclusion 
is continued until < 5.0. Usually, the number of excluded detectors is one 
or a few. In the next step, we search for the core location and the shower size, 
which corresponds to the normalization factor in Equation (7), to maximize 
the likelihood function: 



n 

i=l 



1 



exp 



i=l 



(16) 



where pi is the electron density observed by i-th detector and p{Ri) is the 
particle density estimated from the LDF. The fluctuation of electron density cxj 
takes account of fluctuations in the longitudinal development and the detector 
response. This fluctuation was experimentally expressed by Teshima et al. [25]. 
At this step, we again exclude a detector if its observed density deviates by 
more than 3a, and we assume these signals are possibly overestimated by 
an accidental coincidence or delayed particles. Finally, we estimate the local 
density 5*6) (600) and convert it to the primary energy. 

Figure 13 shows the distributions of energies evaluated using the above anal- 
ysis method for a large number of artificial proton air shower events simulated 
with energies of 3 x lO^^eV and lO^^eV at zenith angles less than 45°. For 
artificial events above lO^^eV and 4 x lO^^eV, 68% have accuracy in arrival 
direction determination better than 2.8° and 1.8°, respectively. These artificial 
events were simulated over a larger area than the AGASA area with directions 
sampled from an isotropic distribution. In this air shower simulation, the fiuc- 
tuation on the longitudinal development of air showers, the resolution of the 
scintillation detectors, and the statistical fiuctuation of observed shower par- 
ticles at each surface detector were taken into account. The primary energy 
is determined with an accuracy of about ±30% at 3 x lO^^eV and ±25% at 
lO^^eV, and the fraction of events with 50%-or-more overestimation in energy 
is only 2.4%. 

Although only events whose cores are located within the array area are used 
in our papers, some events with real cores located near but outside the array 
boundary are reconstructed as "inside" events. The assigned energies of such 
events are smaller than their real values since the core distance of detectors 
become nearer than the true distances. On the other hand, such events that are 
assigned "outside" the array boundary in the analysis procedure against their 
input core locations inside the array are excluded in our selection. The effects 
from these mislocation of cores are taken into account in the distributions in 
Figure 13 and the exposure in Figure 14. 
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4 AGASA energy spectrum and the relation to that in lower en- 
ergy determined at Akeno 



In order to derive the energy spectrum of primary cosmic rays, the obser- 
vation time and the aperture for the selected events must be evaluated as 
a function of the primary energy. The aperture is determined by analyzing 
the artificial showers simulated over an area larger than the AGASA area 
described above. 

Figure 14 shows the energy spectrum observed with AGASA with zenith 
angles smaller than 45° up until July 2002. The exposure (the aperture x the 
observation time) is also drawn in Figure 14 and is almost constant at 5.1 
X lO^^m^ s sr above lO^^eV for the events inside the array boundary. Closed 
circles indicate these "inside" events, and open circles are "well contained" 
events whose cores are located at least 1km inside the array boundary. The 
energy spectra for inside and well contained events agree well with each other 
and hence our criterion of selecting all events inside the boundary can be 
justified. 

Though we have examined the systematic errors in energy determination 
carefully, it is not easy to calibrate the absolute energy experimentally to 
decide whether 10^°eV candidate events really exceed the GZK cutoff en- 
ergy. One method is to compare the spectrum with the extension of Akeno 
energy spectrum measured at lower energies. At Akeno there are arrays of 
various detector-spacing depending on the primary energy of interest, and en- 
ergy spectra have been determined systematically over five decades in energy 
under the similar experimental procedures [29]. 

In the lO^^eV energy region, a comparison of energy determination using 
S'(600) and Ne for each event can be made with the Ikm^ array, where 156 
detectors of Im^ area each are arranged with 120m separation. The relation 
converting Ag to energy at Akeno is determined exprimentally via the lon- 
gitudinal development curve measured at Chacaltaya and Akeno [29] and is 
expressed by 



One of the largest events hitting the Ikm^ array is shown in Figure 15, with 
energies estimated from the shower size Ag and 5'(600). In Figure 16, some 
events are plotted on 5'(600) vs Ag diagram with the 5'(600)-Ae relation from 
Equations (1) and (17). Though the number of events above 10^^ eV is small, 
the difference of energy determined using both methods is within 10%. In 
other words, the energy conversion factor from S'(600) by simulation is in 




(17) 
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good agreement with that from A^e by experiment. 

The energy determined by the Ikm^ array (Ei) and that by the 20km^ array 
(E20), whose detectors are deployed with about 1km separation and is the 
prototype array of AGASA, have also been compared in the lO^^eV energy 
region [30]. The ratio E20/E1 is 1.10 and the dispersion is 45%. Since the 
median energy of the showers is 10^^'^eV, the error in the S'(600) determination 
by the 20km^ array is rather large and hence the wide spread is reasonable. 

In Figure 17, the spectrum obtained by the Ikm^ array (Ei) is shown with 
open squares and that by AGASA by closed squares. There is a difference 
in the overlapping energy region representing a 10% energy difference. In the 
same figure, results below lO^^eV from several experiments are plotted. The 
Akeno energy spectrum is in good agreement with other experiments [31] from 
the knee to the second knee region, except Blanca [32] and DICE [33]. The 
comparison of the present results with other experiments in the highest energy 
region will be made elsewhere. 



5 Conclusion 

We have reevaluated the uncertainties in energy estimation using data accu- 
mulated over ten years. Table 2 summarizes the major systematics and uncer- 
tainties in energy estimation. Here, the symbol "+" means that currently as- 
signed energies should be pushed up under a particular effect, and the symbol 
"— " represents a shift in the opposite direction. The probable overestimation 
of 10% due to shower front structure and delayed particles may be compen- 
sated for by the probable underestimation of the energy conversion factor by 
10%, an effect resulting from the inclusion of the average altitude of AGASA 
and the proper definition of what is meant by a "single particle" . Adding un- 
certainties in quadrature, the systematic uncertainty in energy determination 
in the AGASA experiment is estimated to be ±18% in total. Therefore, the 
currently assigned energies of the AGASA events have an accuracy of ±25% 
in event-reconstruction resolution and ±18% in systematics. 

It should be noted that the Akeno- AGASA spectra cover over five decades 
in energy, connecting smoothly from the knee to a few times lO^^eV, except 
for a 10% difference in energy in the lO^^eV region. This may be due to 
the difference in the energy conversion relations for the experiments and is 
within the systematic errors evaluated here. It is concluded that there are 
surely events above lO^^eV and the energy spectrum extends up to a few times 
10^°eV. The present highest energy event may only be limited by exposure. 
The next generation of experiments with much larger exposures are highly 
anticipated. 
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Fig. 1. An example of the pulse width distribution of a scintillation detector (TB22) 
for one run (about half a day). One channel corresponds to 500ns. This distribution 
is used to monitor the gain and the decay time of the exponential pulse. 
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Fig. 2. An example of time variation of the cable delay (in nanoseconds) between a 
detector (TB22) and the control unit at a branch center. The discontinuity around 
50,000 MJD is due to a change in the detector position and another one is due to 
a system upgrade in 1995. 
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Fig. 3. (a) An example of the time variation in the PWA peak channel (ti) of a scin- 
tillation detector (TB22). (b) Distribution of a{ti)/{ti) for all detectors, determined 
with the data from one month. 
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Fig. 4. (a) Time variation of Aa/a for a detector (TB22). (b) Distribution of 
a{a)/{a) for all detectors, determined with data from one month. 
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Fig. 5. The lateral distribution of energy deposited in an AGASA scintillator (•) in 
units of PH^^,^ is compared with the experimental lateral distribution of AGASA 
(dashed curve). The density of electrons (> IMeV), photons (> IMeV), muons (> 
0.5GeV) and muons (> IGeV) are also plotted. (Left: proton primary; Right: iron 
primary) [Prom Nagano et al., 1998] 
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Fig. 6. Arrival time distribution of charged particles in a 30m^ detector measured by 
a wave form recorder at 1,920m from the core for a 2 x lO^'^eV event. Solid curves 
correspond to to = 800ns. The areas are normalized to the number of particles 
within 2.5/is (87 particles) and 3.5//S (115 particles), respectively. [From Hayashida 
et al, 1994] 
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Fig. 7. Arrival time distribution of shower particles, to for 534m, 750m and 1,050m 
are determined experimentally with a 30m^ scintillator [23], and those for 1,500m 
and 2,000m are extrapolated assuming log to oc r. 
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Fig. 8. Density overestimation due to the effect of shower front thickness estimated 
by considering the LDF profile for different energies. The drop at large core distance 
occurs at a radius where the expected particle count in a detector is one. 
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Fig. 9. Fraction of delayed particles to all shower particles. The current results are 
for to > 3//S and Nij >1.0 with three energy ranges of log(Energy [eV]) = 18.5-19.0 
(big open circles), 19.0-19.4 (big open squares) and above 19.4 (big closed squares). 
The previous Akeno result by Teshima et al. [25] (small open circles), and the result 
by Linsley [26] (small closed circles) are for > 4/is and >3.0 with energies 
around lO^^eV. 
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Fig. 10. Core distance dependence of the ratio Nd/Nz. 
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Fig. 11. Delay time and Nd of delayed particles as a function of core distance. 
Vertical bars indicate the 68% confidence limits in each bin. Data points within 
each bin range from 15 to 30, except for largest core distance bin (3 data points). 
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Fig. 12. Overestimation factor due to delayed particles as a function of core distance. 
The solid line indicates the OF values at each core distance with = 6fis and their 
uncertainties are shown by the shaded region. 
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Fig. 13. Accuracy of energy determination for 3 x lO^^eV and 10^'' eV showers with 
zenith angles less than 45° . The solid curve indicates the Gaussian distribution fitted 
to the histogram. 
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Fig. 14. Energy spectrum determined by AGASA and the exposure with zenith 
angles smaller than 45° up until July 2002. (Open circles: well contained events; 
Closed circles: all events) The vertical axis is multiplied by E^. Error bars repre- 
sent the Poisson upper and lower limits at 68% confidence limit and arrows are 
90% C.L. upper limits. Numbers attached to the points show the number of events 
in each energy bin. The dashed curve represents the spectrum expected for extra- 
galactic sources distributed uniformly in the Universe, taking account of the energy 
determination error. The uncertainty in the exposure is shown by the shaded region. 
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Fig. 15. Comparison of energies determined from and ^(GOO) for one of the 
largest events landing inside the Ikm^ array. 




Fig. 16. ^(GOO) and Ng for the events landing well inside the Ikm^ array. A solid 
line is S'(600)-A'^e relation derived from Equations (1) and (17). 
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Fig. 17. Cosmic ray energy spectrum over a wide energy range. The present AGASA 
energy spectrum is shown by closed squares. The spectrum from the Akeno Ikm^ 
array is shown by open squares. The Akeno-AGASA energy spectrum covers more 
than 5 decades of energy and is in reasonable agreement with most energy spectra 
below lO^^eV. 
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Fig. 14(b). Same plot but only the spectrum. 
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Table 1 



Energy conversion from ^(GOO). The column "Single Particle" describes the defini- 
tion of "a single particle" used in the evaluation of 5(600). Each formula is evaluated 
at the altitude given in the column "Altitude" . 



Simulation 


Single 


Altitude 


Interaction 


Primary 


E = ax 


lO^'^ • 5o(600)^ 


Citation 


Code 


Particle 




Model 


Composition 


a 


h 




COSMOS 


"electrons" 


900m 


QCDJET 


P 


2.03 


1.02 


[15] 


CORSIKA 




900m 


QGSJET98 


P 


2.07 


1.03 


[20] 


(v5.623) 






SIBYLL1.6 


Fe 

P 
Fe 


2.24 
2.30 
2.19 


1.00 
1.03 
1.01 




AIRES 


PKeak 


667m 


QGSJET98 


P 


2.17 


1.03 


[21] 


(v2.2.1) 






SIBYLL1.6 


Fe 

P 
Fe 


2.15 
2.34 
2.24 


1.01 
1.04 
1.02 





Table 2 

Major systematics of AGASA. 
Detector: 

detector absolute gain 

detector linearity 

detector response (box, housing, ...) 
Air shower phenomenology: 

lateral distribution function 

^(GOO) attenuation 

shower front structure 

delayed particles 
Energy estimator 5(600): 

interaction models, chemical compositions (p/Fe), 
simulation codes, height correction, > 
5(600) fluctuation 

Total ±18% 



± 0.7% 

± 7% 
± 5% 

± 7% 
± 5% 

- 5% ± 5% 

- 5% ± 5% 



+10% ±12% 
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